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Abstract 

We study word series and extended word series, classes of formal series for 
the analysis of some dynamical systems and their discretizations. These series are 
similar to but more compact than B-series. They may be composed among them¬ 
selves by means of a simple rule. While word series have appeared before in the 
literature, extended word series are introduced in this paper. We exemplify the use 
of extended word series by studying the reduction to normal form and averaging 
of some perturbed integrable problems. We also provide a detailed analysis of the 
behaviour of splitting numerical methods for those problems. 
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1 Introduction 

In this paper we study word series and extended word series, classes of formal series 
of functions for the analysis of some dynamical systems and their discretizations. We 
exemplify the use of extended word series by studying the reduction to normal form 
of some perturbed integrable problems. We also provide a detailed analysis of the be¬ 
haviour of splitting numerical methods for those problems. Word series are patterned 
after B-series ||23l, a commonly used tool in the study of numerical integrators; while 
B-series are parametrized by rooted trees, word series are parametrized by words built 
from the letters of an alphabet. Series of differential operators parametrized by words 
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(including the Chen-Fliess series) are very common in control theory ll25l and dynam¬ 
ical systems M and have also been used in numerical analysis (see, among others, 
Cll, Ga, Chi). Word series are mathematically equivalent to series of differential op¬ 
erators, but being series of functions, they are handled in a way very similar to the way 
B-series are used by numerical analyst. Word series, as defined here, have appeared 
before in the literature, explicitly d, m or implicitly HI. Extended word series 
are introduced in this paper. 

B-series, introduced by Hairer and Wanner in 1974 ll2^ . provided the first example 
of the application of formal series of functions to the theory of numerical integrators 
(see ll^ for a historical survey). B-series, particularly adapted to Runge-Kutta and 
related methods, give a convenient, systematic way of performing the nontrivial alge¬ 
braic manipulations needed to write the expansion of the local error in powers of the 
stepsize. In addition, they facilitate the construction of integrators found by composing 
simpler integrators; such a construction is required e.g. when investigating the effec¬ 
tive order of Runge-Kutta methods |j71, IS). The usefulness of B-series stems from the 
fact that the composition of two B-series is again a B-series whose coefficients may 
be written down explicitly and are universal in the sense that they are independent of 
the particular differential system being integrated. B-series and their extensions grew 
more important within the notion of geometric integration lf33l . In 1994 Calvo and one 
of the present authors m showed how the conditions for a Runge-Kutta scheme to be 
symplectic may be advantageously derived by examining the corresponding B-series. 
Hairer’s article im started the use of B-series to find explicitly modified systems. 
Since those pioneering contributions the role of B-series and its generalizations EZl in 
geometric integration has kept growing as it may be seen in the treatise ll22l . 

The word series and extended word series considered here are, when applicable, 
more convenient than B-series. A reason for this convenience is that they are more 
compact; in fact the coefficients of a B-series are parametrized by (possibly coloured 
or decorated) rooted trees and there are many more rooted trees with n vertices than 
words with n letters. A second advantage of word series and extended word series over 
B-series is that for words the composition rule (see o and p3| )) is much simpler than 
for rooted trees. 

An overview of the contributions of this paper is as follows. 

Section]^ gives a summary of the rules to manipulate word series. A group Q is 
introduced that plays the role played by the Butcher group in the theory of B-series. 
The solution of the differential system being integrated and some numerical methods, 
including splitting algorithms, may be represented by elements of Q. We also identify 
the Lie algebra g associated with Q and the corresponding bracket. This material is 
very much related to the theory of Hopf algebras ll28l . IQ; however Sectionj^has been 
written with an audience of computational scientist in mind and a number of more 
algebraic considerations have been postponed to Section]^ 

Extended word series are introduced in Sectionj^to cope with perturbed integrable 
problems; roughly speaking we treat problems that may be seen as arbitrary pertur¬ 
bations of systems that, in suitable variables, may be cast in the form {d/dt)y = 0, 
{d/dt)9 = uj (in the language of classical mechanics lH ylO would correspond to 
action/angle variables). We describe the relevant group Q and algebra g. 

In Section|^we show how to use extended word series to bring perturbed integrable 
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problems to normal form, i.e. how to change variables to reduce the system being 
analyzed to a form as simple as possible. As distinct from standard ways of finding 
normal forms, the extended word series approach does not rely on the vector field 
being polynomial. Furthermore, the computations required here are universal (in the 
sense of HI): they are independent of the particular system under consideration. 

For highly oscillatory problems the reduction to normal form is very much related 
to the process of averaging out the oscillatory components of the solution and therefore 
Sectionj^extends the material in the series of papers ifT^ . lITSll . 1141 . iflSl . Furthermore 
normal forms readily lead to the explicit computation of (formal) invariant quantities 
of dynamical systems and their discretizations, an issue not covered in this paper and 
treated in the follow-up article ll30l . 

Section]^ is devoted to the study of general splitting algorithms to simulate per¬ 
turbed integrable problems. We show how our algebraic approach leads to a convenient 
expansion of the local error. It is well known that the behaviour of the corresponding 
global error as h varies is unfortunately extremely complex, as it depends on arithmetic 
relations between h and the periods present in the solution. Extended word series pro¬ 
vide a powerful instrument to analyze that behaviour. In fact, two different approaches 
are put forward here. In the first, the integrator is processed, i.e. subjected to changes 
of variables, to remove oscillatory components. The second approach is based on con¬ 
structing a modified system for the integrator and then bringing the modified system to 
normal form. Of much interest is the fact that the validity of the modified system holds 
even if h is not small relative to the periods in the problem (cf the use by Flairer and 
Lubich of modulated Fourier expansions jlSSl l. The techniques in this section may be 
readily applied to the construction and analysis of improved integrators, such as those 
considered in e.g. ED, 01; this will be the subject of future work. 

Sectionj^contains proofs and technical material and there is an Appendix devoted 
to the practical applicability of splitting integrators. 


2 Word series 


This section presents word series and provides a summary of the rules for their appli¬ 
cation. The presentation has computational scientist in mind and focuses on essential 
features; additional details and proofs are given in Section |6T| where the approach is 
more algebraic. Until further notice all functions are assumed to be smooth. 

2.1 Definition of word series 

We consider the £)-dimensional initial-value problem given by 


a;(0) = xo 


( 1 ) 


and 



( 2 ) 
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where t is the (real) independent variable, A is a finite or infinite countable set of 
indices and for each a G A, Xa is a scalar-valued function and fa a I?-vector-valued 
map. 

It is well known that the solutions of Q may be expanded formally as follows. As¬ 
sociated with each vector field fa in Q, there is a first-order linear differential operator 
Ea'. if p is a scalar-valued function, then the function Eag is defined by 

D „ 

Kgix) = XI (x) g^9{x) (3) 

f=i 

(superscripts denote components of vectors). We shall also let Ea act on vector-valued 
mappings; it is then understood that the operator is applied componentwise. If x{t) 
satisfies the chain rule yields 

^g{x{t)) = X Xa{t){Eag){x{t)) 


or 


g{x(t)) = g{x(0))+ ^ / dtiX{h){Eag){x{h)). 

The same procedure maybe now applied with {Eag){x{ti)) in lieu of g{x{t)) to rewrite 
the last equation as 


9{x{t)) 


9{x{0)) 


aeA ‘ 


dti X{ti){Eag){x{ 0 )) 


+Exr dtiXiti) f dt2\(t2){Eb{Eag))[x{t2)). 

aeAbeA-^^ -^0 


By continuing this Picard iteration and setting g equal to the identity function g{x) = x, 
we find that the solution of Q-Q has the formal expansion 

OO 

x{t) = xo + X X (f )/*ai ■■-a,i(xo) , (4) 

n—1 ai,...,a,iGA 


where the vector-valued mappings fa^...a„{x) and the scalar-valued functions am--a, 
satisfy the recursions 


fai---an{x) — dxfa2---an{x) fafxf TT > 1, 

(dxfa 2 ---a„{x) denotes the value at x of the Jacobian matrix of /as -a^) and 

*Aai(f) — / Xai{tl')dti^ 

Jo 

— / Xa^ifn) (y-ai---an — Jfn) dtyi^ 72 > 1. 

Jo 


(5) 

(6) 
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For future reference we note that 




^ani^n) 




rt2 

i(^n—l) ' I dt\ 
Jo 


or 



Aai(fl) ’ ’ ‘ dti ' ' • dtj^j 


where the n-fold integral is taken over the simplex 
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Sn{t) = {(fi,..., in) G K” : 0 < fi < • • • < < t}. 

Let us present some examples (more may be seen in 

1. In the simplest illustration, the set A has only one element a and the corre¬ 
sponding A takes the value 1 for each t. Then ([2| is the autonomous system 
{d/dt)x = fa(x)- For each n, the inner sum in ([^comprises a single term and 
from the corresponding coefficient is found to be aa---a{t) = t^jn\. In this 
case (OT is the standard Taylor expansion of x(t). 

2. The autonomous system {d/dt)x = F{x) with the right-hand side split as F{x) 
= fa{x) + fb{x) is of the form ([^ with A = {a,b} and Xa{t) = Ab(t) = 1. 
For each n, the inner sum in Q comprises 2” terms and each of them has a 
coefficient f"/n!. The expansion is the Taylor series for x{t) written in terms 
of the pieces fa and /& rather than in terms of F, a format that is useful in the 
analysis of splitting numerical integrators. It is of course possible to split F in 
m > 2 parts or even in infinitely many parts; then A has m or inhnitely many 
elements. In all cases the integral (j^ has the value jn]. 

3. Let w > 0 be a fixed number. If A = Z (the set of all integers) and, for k G 
Z, Afe(f) = exp(zkcot), then (|^ is a non-autonomous system 27r/a;-periodic in 
t with the right-hand side expanded in Fourier series. Here the functions Xk 
possess complex values; if the system Q is real then the mappings take values 
in and fk is the complex conjugate of f-k for each k G Z. The expansion 0 
has been used in ifT^ . lfT3l . ifT^ . ifTSl to study analytically periodic problems. 

4. Also treated in OS, 111,1113 are quasiperiodic problems. If w G is a vector 
of frequencies, these are of the form (j^ with A = Z‘^ and 

Ak(f) = exp(ik • wf), k s (8) 


This case will be taken up in the next section. 

The notation in Q may be made slightly more compact by considering A as an 
alphabet and the strings oi • • • a„ as words with n letters. Then, if yV„ represents the 
set of all words with n letters, 0 reads 

OO 

x{t) = Xo + '^ ^ a^u{t) fw{xo)- 

n—1 TtlGWn 
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If we furthermore introduce the empty word 0 and set Wg = {0}, ait = 1, = 

X, then the last expansion becomes 

OO 

= ^ aw{t)fJ^XQ)= ajt) fj^xo), (9) 

n—0w£Wn wGW 

where W represents the set of all words. 

Subsequent developments will make much use of the set of all mappings i5 : 
yy — > C; if 5 S and w is a word, then 6^^ is a complex number. This set is 
obviously a vector space for the usual operations between maps: if /ii, /i 2 are scalars 
and(5i, ^2 G C^, then (/ri(5i+/r2<52) G is defined by {fiiSi+fi 2 S 2 )w = Mi(<^i)i« + 
tJ-2{S2)w for each w eW. 

The expansion (j9]) motivates the following definition: 

Definition 1 If 6 G C^, then its corresponding word series is the formal series 

Ws{x) = Y, Sn,Ux). (10) 

The scalars and the functions fyj will be called the coefficients of the series and 
word-basis functions respectively. 

Clearly the word-basis functions change with the mappings fa in the system (|^ 
being studied. With this terminology, for each fixed t, the formal series 0 for the 
solution x{t) of is a word series whose coefficients aw (t) are given by 0 and 

a^{t) = 1 (these coefficients are independent of the mappings fa). Also, for each t, 
the right-hand side of 0 is a word series with coefficients /3a{t) = Xa(t) for words 
with one letter and / 3 u,(f) = 0 for all other words. As we shall see later, word series Ws 
corresponding to other choices of coefficients are useful in the analysis of dynamical 
systems and their numerical integrators. 

Remark 1 Word series and moulds. In Ecalle’s terminology, word series are moulds, 
see im, Ga, and the convolution product considered below is the mould product. 
(Ecalle ini also considered a composition of moulds —not discussed in this paper— 
that is analogous to the substitution of B-series ini.) 

Remark 2 Word series and B-series. Each fw, w built up from partial deriva¬ 

tives of the fa, a G A, e.g., if a,b,c G A, 

fba{x) = da:fa{x) fb{x), 

fcba{x) = dxfba{x)fc{x) = d^x fa{x)[fb{x), f^x)] -f dxfa{x) dxfb{x) fc{x). 

The functions dxfa{x) fb{x), dxxfa{x)\fb{x),fc{x)],dxfa{.x) dxfb{x) fdx) in these 
expressions are examples of elementary differentials; each word basis function fw{x), 
with w G Wn, n > 0, is a linear combination with integer coefficients of elementary 
differentials of order n (i.e. containing n functions fa{x)). There is an elementary 
differential corresponding to each A-coloured (or A-decorated) rooted tree, i.e. to each 
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rooted tree where to each vertex it has been assigned an element of A. By expanding 
each word basis function in terms of elementary differentials, the series ( [T0| ) becomes 
a so-called B-series 

T 

where the summation is extended to all A-coloured rooted trees and (a;) is the ele¬ 
mentary differential corresponding to r. B-series were introduced by Hairer and Wan¬ 
ner ll2^ in the simplest case where the alphabet A has only one letter. B-series corre¬ 
sponding to this and larger alphabets are often used in numerical analysis; word series 
being more compact are better suited to analyze some integrators. For the relation be¬ 
tween word series and B-series see ll28l and lITSl (Ecalle used in this connection the 
terminology arborifaction-coarborification). 


Remarks Word series as power series. For a system {d/dt)x = ^a{t)fa{x), 

where e is a scalar parameter, becomes the formal power series 

OO 

n—0 wGWji 


Note that when the alphabet A is infinite the coefficient of e” is itself an infinite series 
that has to be understood formally. The format (10 1 is of course recovered from the 
power series by setting e = 1; therefore both formats are equivalent. While the papers 
113, 03 use the power series format, we prefer to work with ( [T0| ) as it leads to more 
compact formulae. Some readers may find it useful to mentally substitute efa for 
fa everywhere in what follows; this may be particularly the case for the perturbed 
integrable problems considered in Section |3 


2.2 Operations with word series 

2.2.1 The convolution product 

Given 6, 6' G C^, we associate with them its convolution product S-kS' G defined 
by 

n—1 

i=i 

(here it is understood that (5 * (5')0 = The convolution product is not commu¬ 

tative, but it is associative and has a unit (the element 1 G with 10 = 1 and 
1 u) = 0 for w 7 ^ 0). 

As we shall see, the operation ★ plays an essential role in the manipulation of word 
series. 


2.2.2 The group 

If w G yVm and w' G Wn are words, m,n > 1, their shuffle product wluw' EH 
is the formal sum of the (m + n)!/(m!n!) words with m + n letters that may be 
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obtained by interleaving the letters of w and w' while preserving the order in which 
the letters appear in each word. (Examples: for the words ab, cd, the shuffle product is 
abujcd = abed + acbd + cabd + aedb + cadb + edab, for the words ab, a, the product 
is a 6 Lua = aba + aab + aab) In addition 0Ljjr(; = r(;Ljj0 = uifor each w G W. The 
operation m is commutative and associative and has word 0 as a unit. 

We denote by ^ the set of those 7 G that satisfy the so-called shuffle relations: 
70 = 1 and, for each w, w' G W, 


^w^w' — 


N 

i=i 


if 


N 

i=i 


( 12 ) 


The set Q with the operation * may be regarded in a formal sense (cf. 0 ) as a non- 
commutative Lie group (see Section 6.1.2| l. For each fixed t, the family of coefficients 
defined by (j^ and a$ (t) = 1 is an element of the group Q (to prove this, consider ( 12 1 
for w G Wm and w' G W„ and use Q to write ayjCtw' and each aw^ as integrals over 
subsets of cf. EU Corollary T5]). 

When 7 belongs to Q, the word series {x) has properties that are not shared by 
general word series. For j G G, changes of variables x = C{X) commute with the 
formation of word series as described in ifTSl Proposition 3.1]. Moreover, for j G G, 
W„f{x) may be substituted in an arbitrary word series Ws{x), 5 G to get a new 
word series; more precisely 


Ws(W^(x)) = Wy^s(x), 


(13) 


i.e. the coefficients of the word series resulting from the substitution are given by the 
convolution product 7*(5 (this is proved in Section 6.1.3[ ). A similar rule exists of course 
for B-series, but the recipe there is more complicated than (111 ||23l, ll22] Chapter III]. 

In the numerical analysis of differential equations, word series with coefficients 
in G appear e.g. as expansions in powers of the stepsize of splitting integrators, see 
Section Then 0 provides the recipe to compose integrators or to compose an 
integrator and a mapping. For each fixed t, the right-hand side of (j^ is an example of 
a word-series with coefficients in the Lie algebra p of the group G that we study next. 


2.2.3 The Lie algebra g 

We denote g the set of elements /3 G such that /3^ = 0 and for each pair of 
nonempty words w, w', 

N N 

Pwj = 0 if wLuw' = Wj. 

i=i i=i 

It is clear that g is a vector subspace of the vector space C^; furthermore g is closed 
for the skew-symmetric product defined by 

[/?,/?'] ( 14 ) 
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This product satisfies the Jacobi identity and therefore endows g with a structure of Lie 
algebra (see Section 6 .1.2 1 . In fact g is the Lie algebra of the Lie group Q: the elements 
/3 € g coincide with the velocities at 1 G G of curves in G- In symbols, if 7 (f), t G M. 
is a curve in G such that 7 ( 0 ) = 1, then /? G defined by 


/3 



t=o 


(15) 


(i.e. /3u, = (d/dt)'yiu(0) for for each w G W) belongs to g. Moreover any /3 G g arises 
in this way: the exponential 


exp*(f/3) = 1 +51 7!^*^ 

J=i 


(16) 


(/3*-^ is the convolution product of j factors all equal to /3) defines a curve of elements 
of G with velocity /3 at f = 0. The points of this curve actually form a one-parameter 
subgroup of G since exp,^(f/3) *exp^(f'/3) = exp^((t + The exponent /3 may be 
retrieved from the exponential 7 = exp^(/3) by means of the logarithm 


P = log* ( 7 ) 




1 )*+ 


(17) 


Just as the convolution product with an element of G corresponds to the opera¬ 
tion of substitution of the associated word series (see ([T3]l), the convolution bracket 
( [T4 ) i corresponds to the Jacobi bracket (commutator) of the associated word series, for 

P,P' G g: 


{d^Wij>{x))Wij{x) - {d^Wi3{x))Wi3i{x) = 


To prove this, let 7 (f) be a curve with velocity /3 as above; then, for any 5 G 


(da:lVg(x))Wf3(x) 


^W4w,(,)(x)) 


d 

dt 






W0^six). (18) 


(We have successively used the chain rule, 0, and the bilinearity of *.) 

Since for P G g, the word series W/ 3 (x) belongs to the Lie algebra (for the Jacobi 
bracket) generated by the mappings (vector fields) fa, the Dynkin-Specht-Wever for¬ 
mula II 24 I may be used to rewrite the word series in terms of iterated commutators of 
these mappings: 


LXJ - 

= 51 /3ai-a„[[---[[/ai,/a2],/a3]---],/a„](a;). (19) 

n—1 

(For n = 1 the terms in the inner sum are of the form Pa^ fai (x)-) 
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2.2.4 Nonautonomous differential equations in Q 

Initial value problems 


—x{t) = a:( 0 ) = xo, 


dt 


( 20 ) 


where for each t, f3{t) G g, are a natural generalization of (we recall that the 

right-hand side Q does not include contributions from basis functions associated with 
words with more than one letter). These problems may be solved formally by using the 


ansatz x{t) = Wa(t)ixo) with a(t) G Q for each t. In view of (13 i, we may write 


dt 


^^a(i)(^o) (^^a(t) (^o)) (^o)i ^^a(0)(^o) ^0: 


which leads to the linear, nonautonomous initial value problem 


— a(f) = a(f) 


a(0) = a. 


( 21 ) 


For the empty word, according to (11 1 , {d/dt)a% = a 0 (f)/ 30 (t); since /3(f) G g implies 


/ 30 (f) = 0, we see that a 0 (f) = 1. For a word a G Wi, {d/dt)aa{t) = a 0 (f)/ 3 a(f) -f 
aa{t)P(/t{t), which leads to aa{t) = /q/ 3a(fi) dfi- The process may be continued 
in an obvious way and induction on the number of letters shows that 0 uniquely 
determines aw{t) for each w G W. Furthermore, for each f, the element Q!(f) G 
defined in this way belongs to Q-, while this may be established by means of the Magnus 
expansion (see e.g. Q, ll22l Chapter IV]) we provide an elementary proof in Section]^ 
Conversely, any curve a{t) of group elements with a(0) = 1 solves a problem of 
the form ( |2T| with 

' d 


/3(f) = a(f) 


-1 




Since 


>i{ty 


dt 


^*a(f+ s)) 


s=0 


for each f, the element /3(f) defined in this way is a member of g. 

The investigation of normal forms below is based on changing variables. A change 


of variables x = IFk(-A), k G G, transforms the problem (20 1 into 




2f(0) = Xo, 


with B{t) -k K = K-k P(t) (or B{t) = Kk /3{t) k K ^), Xo = W ^-1 (xo) (k ^ is the 
inverse of k in the group G)', this is a direct consequence of ([T3]l and (fTS]). 


2.2.5 The Hamiltonian case 

Consider now the particular case where the dimension of (j^ is even and each fa{x) 
is a Hamiltonian vector field llTSl . i.e. fa{x) = J~^\/Ha{x), where J~^ is the standard 
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symplectic matrix. Recall ||2l that the Jacobi bracket (commutator) A, 

of two Hamiltonian vectors fields is again a Hamiltonian vector held and that the cor¬ 
responding Hamiltonian function is the Poisson bracket of the Hamiltonians A and B, 
dehned by {A,B}(x) = VA{x)'^J~^B{x). According to (19i, for each /3 G g, the 
vector held Wp{x) is Hamiltonian 

Wp{x) = j-^vnpix) 


with Hamiltonian function 


H/six) = ^ P-ujH^ix), 

luG W, 

where, for each nonempty word w = ai ■ ■ ■ an, 

Hu,{x) = ■ ■ ■ },HnJ{x). ( 22 ) 

For Hamiltonian systems, changes of variables x = Wk,{X), k G (/, are canon¬ 
ically symplectic; after the change of variables the system is again Hamiltonian and 
the new Hamiltonian function is obtained by changing variables in the old Hamiltonian 
function El. 


3 Extended word series 

In this section we adapt the preceding material to cover perturbed integrable problems. 


3.1 Perturbed integrable problems 


We now consider systems of the form 


y 


'O' 

9 


UJ 




where y G 0 < d < ZJ, w G is a vector of frequencies tvj > 0, j = 1,d, 

and 9 comprises d angles, so that f{y, 9) is 27r-periodic in each component of 9 with 
Fourier expansion 

/(d: ^) = X! /k(d) 

kGZ'i 

ifk{y) and f-ii(y) are mutually conjugate, so as to have a real problem). Systems 
of this form appear in many applications, perhaps after a change of variables (see the 
Appendix). When / = 0 the system is integrable (the angles rotate with uniform 
angular velocity and y remains constant) and accordingly we refer to problems of this 
class as perturbed integrable problems and to / as the perturbation (some readers may 
prefer to substitute ef for /, see Remarkj^. 

After introducing the functions 

hiy,0) =exp(i'k-9) fk{y), y G 6> G (23) 
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that satisfy the fundamental identity 


fk{y, Oi + 02) = exp(ik • 01 ) /k(y, 62), 
the system takes the form 

kez<i 

To find the solution with initial conditions 


dt 


fiy,o) = 


2/(0) = 2/0, ^*(0) = 00, 


(24) 


(25) 


(26) 


we perform the time-dependent change of variables 0 = 77 + fw to get 


^ \y 

dt t] 


^ exp(ik-a;i)/k( 2 /,? 7 ), 

kGZ'' 


(27) 


a particular instance of the problem considered in Section]^ The alphabet A coincides 
with Z'^, and, for each Tetter’ k, Ak(f) is given by The formula (j^ yields 


'y{t) 


y(o)' 

_rjit)_ 


_ 0 (O)_ 


+ E E «ki...k„ (o/ki-k„ (2/(0), 0(0)), 

n—l ki,...,kTi 


(28) 


where the coefficients a are still given by Q (but recall that now the letters a are 
multiindices k) and the word basis functions are defined by (j2^ and Q (the Jacobian 
in (j^ is taken with respect to the U-dimensional variable {y, 0)). We conclude that, in 
the original variables, the solution flow of ([25]l, has the formal expansion 


(t>t{yo,yo) 


y(i)' 


2/0 

_ 1 _ 

'o' 



.^0. 




+E E «k,..k„ it)fk,-kAyo,0o)- (29) 

n—l ki,...,k„ 


Note that the word basis functions are independent of the frequencies lo and the coeffi¬ 
cients a are independent of f. Also from ( [24) l we have the identity; 

fki-k„{y,0i + 02) = exp(i(ki H-h k„) • 0i) /ki-k„ (2/, (* 2 ), (30) 


and, in particular 


/ki -k„( 2 /,^) = exp(7(ki -f-[-k„) • 0) /ki...k„(2/,0), (31) 


With the notation of Section]^ we write 

a: = (2/,0)): 

x{f) = ^ ® 


tuj 


in the following form (here and later 

Wo,{t){xo). 


In order to make the formula even more compact, we introduce the vector space C = 
£d 0 £W (Jefjjje: 
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Definition 2 If (v,S) € C, then its corresponding extended word series is the formal 
series 


(v,S) {^) 




f ■ 


Then the solution ( [29] l of ([25|)-(p6|) has the expansion 

x{t) = W(^t^^ait)){xo), {t 0 J,a{t)) e C, 

with a{t) G G C as defined in Section]^ Also the right-hand side of (25 i is an 
extended word series with /? € g C defined by 

/3^ = 1 if w e Wi, /3„ = 0 if w ^ Wi- (32) 


3.2 Operations with extended word series 

3.2.1 The operation ★ 


The following two linear operators will appear repeatedly. If u is a d-vector, is the 
linear operator in that maps each S G into the element of defined by 

{Ey6)$ = 6$ and 

(E;„d)u, = exp(j(ki H-|-k„)-u)du,. (33) 

for u> = ki... k„. The linear operator on is defined as follows: = 0, 

and for each word u> = ki • • • k„, 

{^vS)w = i(ki H-h k„) • u (5i„. 


Thus Ey and are diagonal operators with eigenvalues exp(i(ki + • • • + k„) • v) and 
i(ki + • • • + k„) • V respectively. Observe that S„(7*i5) = (S„7)*(S„(5) if 7 , d € 
and that: 


d„ 
dt * 


^tv^v Cv 


The symbol G denotes the subset of C comprising the elements (u, 7 ) with u G C‘^ 
and 'j G G- For each t, the solution coefficients (fw, a{t)) G C found above provide an 
example of element of G- With the help of we define an operation as follows. If 

(u, j) G G and {v, 6) G C, then 


(u, 7 )^(w, S) = ( 70 U + Sun, 7 * (S^d)) e C. 


By using (13 1 and (30 1 , it is a simple exercise to check that G acts by substitution on 
extended word series as follows: 


W^0.5)(W^(«,7)(a^)) = jGG- (34) 

In fact we have defined the operation ★ so as to ensure this property. The set is a 
group for the product -k and and G may be viewed as subgroups of The unit of 
G is the element 1 = ( 0 , 1 ). 

* Consider the group homomorphism from the additive group C'* to the group of automorphisms of Q 
that maps each /r S into . Then Q is the (outer) semidirect product of Q and the additive group 
with respect to this homomorphism. 
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3.2.2 The Lie algebra g 


As a set, the Lie algebra g of the group Q consists of the elements {v, S) G C with 
5 G p. Let us describe the bracket in g. Given v G and <5 G g, we have the trivial 
decomposition 


^{v,5) = 


+ Ws (x) — kL(u,o) + ^(o.<5) ■ 


By using (31 1 , one can check that the Jacobi bracket of the vector helds kL(A,o) 
W^(o,5) is 


[W^(w,0), W^(0.5)] = 

From these relations we conclude that, for arbitrary {v, S), {u, p) G 0 g, the Jacobi 
bracket of the vector helds is given by 


=W ^(0 

Accordingly the bracket of g has the expression 

The 0 rehects the fact that is an Abelian subgroup of Q. 


3.2.3 Nonautonomous differential equations in Q 

The initial value problem 


—a:(f) = VL(^,^(t))(a:(f)), x(0) = xq, 

where {uj,i3{t)) G g for each t, may be formally solved in a manner that is exactly 
parallel to treatment given above to (20i: x{t) = W(^tui^a{t)){xo), where q;(0 ) = 1 
and 

Observe that the right-hand side of this equation is, by dehnition of ★, equal to 
(fw, a{t) -k (St^/3(f))), so that a{t) is the solution of an initial value problem of the 
form (21 1 with /3(f) replaced by 

A change of variables x = K)(Ar), (u, k) G Q, can be seen to transform the 

problem into 

jX{t) = ^(0) = Xo, 

where now -B(f) is determined from 


B(t) * K + = K * (^„^(f)) 

and, of course, xq = or Xq = kL(„^„)-i(xo). (Note that {v,k)~^ = 
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3.2.4 Perturbed Hamiltonian problems 


To end this section, assume in (25 1 , that the dimension D is even with D/2 — d = m > 
0 and that the vector of unknowns takes the form 


X = {y, 0) = (p\ ... g™; a\ ..., ( 


where is the momentum canonically conjugate to the co-ordinate and is the 
momentum (action) canonically conjugate to the co-ordinate (angle) 9^. If each /k(a:) 
in is a Hamiltonian vector field with Hamiltonian function then the system 

(|25|l is itself Hamiltonian for the Hamiltonian function 


d 

j = l kGZ'^ 


For each (w, /3) G 0 , the extended word series is a Hamiltonian formal 

vector field, with Hamiltonian function 






(35) 


with Hy^(x) as in (22i. Note that the Lie bracket in g can be used to compute the 
Poisson bracket of formal Hamiltonian functions of the form i 


4 Normal forms and averaging 


In this section we show how the algebraic machinery introduced above may be ap¬ 
plied to build a theory of normal forms m,i3a for the perturbed integrable problems 
of the form (p5[l. This theory hinges on the fact that the linear operator W (o, 5 ) 


\W (c^,o)) W ( 0 , 5 )] (<5 G g) coincides, as we have seen in Section 
operator IL(o, 5 ) i-G IL(o,{,^ 5 )- 

Let us consider an autonomous system 


3.2.2 


with the diagonal 


dt^ dt 


y 


'o' 

_e_ 


UJ 


+ Wfi{x) = W(^^fi){x), /3 G g. 


(36) 


As noted before, this format yields the perturbed problem (25 i when (3 is chosen as 
in ( [32| ). The more general case where /3 is any element in g will be necessary to deal 
with splitting integrators later. We shall change variables x = W,^{X) = IL(o,k)(-V), 
tv € (I, in order to simplify ([3^ as much as possible. 


Remark 4 There is nothing lost by assuming that x = Wk{X) is a (not extended) 
word series in the new variables X —or equivalently an extended word series of the 
special format x = W (o,k) (3f)—■ More general changes x = W („ ^) (X) do not allow 
for additional simplications in (j3^. 
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From Section 3.2.3 we know that the transformed system is 


—X=- 
dt dt 


y 


'o' 

0 


UJ 




with 


+ j5 -k K = K-k p. 


( 37 ) 


(38) 


Our aim is to choose fd G q and k G G subject to ( 38 1 and such that (3 is as simple 
as possible; then the system is said to have been brought to normal form. Of course the 
maximum simplification would be obtained by setting ,0 = 0, but for this choice of fd 
it is not possible to find an appropriate k; this will be clear in the proof of Theorem 
and is to be expected from general results on normal forms 111 , 1321 . More precisely, 
perturbations that commute with cannot be eliminated by changing variables. 

One then has to restrict the attention to /? € g such that in (371 the unperturbed vector 
field and the perturbation commute, i.e. [kF(a;p), kF^g = 0. This is equivalent to 
kF(g j ^) = 0, or, in terms of the coefficients, 

z(ki + • • • + k„) • (jj /3ki...k„ = 0, (39) 

for each nonempty word ki... k„. We have then the following result, which is proved 


constructively in Section 6.2 


Theore m 1 There is a change of variables x = Wf^{X), n G G, that reduces the 
system (36) to the form (37), where Id G Q and = 0 for all words w = ki ... k„ 
such that (ki + • • • + k„) ■ w f 0. Furthermore the vector fields W (^o){X) and 
kF|.g commute and the solutions of (37) satisfy 


X{t) = fit(^X{0) 


+ 


= Mxm 


where fit is the solution flow of the system (d/df)X = Wp{X). Equivalently, X{t) = 
W(tc,,ait)){XiO)), where 

{tuj,a{t)) = {toj,expfytfi)) = expfytfi)^{tuj, 1) = (fw, l)^exp^(f^). 


If the system ( [36| ) is Hamiltonian, the change of variables is canonical symplectic 
and the transformed system is Hamiltonian. 

When w is nonresonant, i.e. k • w 7^ 0 for k 7^ 0, the theorem implies, in view 
, that the transformed vector field independent of the angular 

es 0. In other words, (371 is a system where the angles have been averaged HI, 


of (31 
variaF 

Cl, fyiX . In the general situation with a nontrivial resonant module 

= {k e : k • w = 0}, 


the transformed vector field depends on 0. However this dependence is only through 
a number of combinations li • 0 ,..., 1^ • 0, r < d, where li, ..., 1^ Gif are linearly 
independent and span the resonant module. 
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Remark 5 Consider the highly oscillatory case where ( [36] l depends on a small param¬ 
eter S and Lo = 0{l/5), Wp{x) — 0{1). The combinations not eliminated by the 
change of variables have the property that their velocities {d/dt)li ■ 0 are 0{1), as dis¬ 
tinct from the situation for the original angles with {d/dt)0 = 0{l/5). In this sense, 
the/flsf angles have been averaged when forming ([T7|. 


For convenience, we shall use the expression oscillatory word to refer to those 
words ki... k„ for which (ki -f • • • -f k„) • cu ^ 0. Thus the theorem may be rephrased 
as saying that the contributions to the vector field corresponding to oscillatory words 
may be removed from ([36| by means of a change of variables. Note that the set go 
of all /3 S g such that /3u, = 0 for all oscillatory words is a Lie subalgebra of g; this 
follows from the fact that /3 is in go if and only if (w, 0) and (0, /3) commute. 

If we now express the commuting vector fields o)(-^) ^)(^) terms 

of the original variables x by applying the recipe for changing variables given in Sub- 
section [J.2.3 we obtain a decomposition of the right-hand side of (361 as a commuting 
sum of two terms: 


The second of these generates a flow 

X{t) = 

where the motion of the fast angles has been averaged. The former generates a quasi- 
periodic flow 

Finally, the flow of ( [36] l is given by 

Remark 6 In the nonresonant case, the commuting decomposition of W has 
been obtained in ifTJl Theorem 5.5] by means of a different (but related) technique. 


5 Splitting methods 


Splitting algorithms llTSl . Il22ll are natural candidates to integrate perturbed integrable 
problems. In this connection, it is extremely important to emphasize that the practical 
implementation of splitting methods is not necessarily based on the simple format ( [25] l. 
Such a simple format is typically reached after suitable changes of variables and is quite 
convenient for the analysis. These points are discussed in the Appendix. 

Given real coefficients, Qj and bj, j = 1,... ,r, we study the splitting integrator 


for (25 I defined by 


4>h = 


du) 

rh ' 


o 6^ 'O' 
't'a ^ 


AP) AU) 

°Kh°Ah- 


(40) 
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Here h is the step-length, the mapping in that advances the numerical solution 
over one time step, and and denote respectively the exact i-flows of the split 

systems corresponding to the unperturbed dynamics 


y 


'O' 

e 


UJ 


(41) 


and the perturbation 


If we set 


d \y 

Jt e 


fiy,o). 


= E' 




3 = 1 


(42) 


the integrator is consistent if a = b = 1. 

Since the unperturbed dynamics with frequencies ujj is reproduced exactly by ( |40) l, 
one would naively hope that the accuracy of the integrator would be dictated for the 
size of / uniformly in oj. It is well known that such an expectation is unjustihed, see 

e.g. m, M- 


5.1 Extended word series expansion of the local error 


Clearly, the mapping has an expansion in extended word series 

= W(^tub,i)ix), (fw, 1)6^; 

furthermore, using Example 2 in Section]^ 

(0,T(f)) G Q, 

where r(f) G Q comprises the Taylor coefficients, i.e. Tu,(f) = f"/n! if tu G Wn- 
The following result makes use of the algebraic formalism to provide explicitly the 
expansion of the numerical solution. 


Theorem 2 The splitting integrator (jj^ in ( 40 • possesses the expansion 

4^h{x) [hauj ^a{h))ix') t 

where aih) G Q is specified by a${h) = 1 and, for n = 1,2,..., 


aki-k„(^) = /i” 
Here, 
and, 

o-i. 


E 




Cj — Ui + • • • + Oj, 
1 


exp(i (cjiki + • • •-f Cj^k„) • wh). (43) 


1 < j < r, 


d 


if 


_ 1 

— ^ <Xje+i---3„ 


if 


<n, jl = ■■■ = jl < ji+l <■■■ <jn- 
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Proof: From (341, we know that (f>h has an expansion in extended word series and that 
the family of coefficients is given by (pay attention to the ordering) 


(oi/iw, l)'k{0,T{bih))'k ■ ■ ■ 'k{arhuj, TL)'k(0,T{brh)); 

it is enough to compute, according to the definition, the products if in this expression. 

□ 

Remark 7 The associated quadrature rule. For words with one letter, the theorem 
yields: 

5k(^) = h bj exp(icjk • ujh). 

l<_ 7 <r 

This obviously corresponds to the approximation of the exact coefficient ak(^) (de¬ 
fined in •S' and (|^) by the (univariate) quadrature rule that on the unit interval has 
abscissas Cj and weights bj. This rule will be consistent if 6 = 1, which is implied by 
the consistency of the integrator. 

Remark 8 Associated cubature rules. Similarly, for n > 1, ( |4^ corresponds to ap¬ 
proximating Q with a cubature rule for the simplex. If the univariate quadrature is 
consistent, so is the cubature rule for each n > 1, because, by using the multinomial 
expansion. 


E 


• • • 5^’- 


• • • ^3^ ^ _ 

cr,. ...i ^ ni! • • • rir-! 

niH- \-nr—n 


= )^(E^^E=(E^^)"Vo1(5(1)). 


i=i 


i=i 


Discussions perhaps become clearer by introducing scaled coefficients and 
such that 


aki - k„(^) = /i"^ki-- k„(^), Q:ki - k„(^) = ^”^ki -k„ (M- 
Note that, by performing the change of variables tj = htj, j = 1,. .., n, in (0, 

^ki-k„(/i) = / ■ • • / exp(i(ykiH- \-t'^kn) ■ u}h) dt[ ■ ■ ■ dt'^, (44) 

J Js„(i) 


and that, therefore. 


|< Vol(5„(l)) = 


With these preparations, we have proved our next result: 


Theorem 3 The local error of the splitting integrator fh in {40) possesses the expan¬ 
sion 

fh{A) fhitc) — VF [h{a—l)u)Jx{h) — OL{h)){A). 
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i.e. 


4>h{xQ) - (j)h{xo) = 


h{a — l)t 


(45) 


+E'*” E (^ki---kTi {^) ^ki-'-k^i (^)) /ki---kTi (^o)■ 

"=1 ki,...,k„GZ‘* 


5.2 Estimates 

In order to obtain error estimates it is now necessary to truncate the infinite series in 


and we shall do so in the next theorem, whose proof is given in Section 6.3 We 
assume hereafter that; 

1. The function f{x) = f{y, 6) is defined in a set = Bfi{yo) xT'^, where Bji{yQ) 
is the ball {y '■ \y — yo\ < R} C 


2. There exists a finite set of indices X C such that for k 
coefficient vanishes. 


X the Fourier 


3. There exists an integer N > 2, such that the Fourier coefficients and their 
partial derivatives of order < — 1 are continuous and bounded in Bii{yo). 

In the first of these hypotheses the form of the domain 12 is natural since f{y, 9) is 
periodic in each of the components of 9. As shown in e.g. Qa, the second hypothesis 
may be relaxed for the conclusions of the theorem to hold; it makes however possible to 
avoid distracting technicalities. In this connection it should be noted that, for nonlinear 
problems, even if / has a finite number of Fourier modes the solution x{t) in (291 will 
include arbitrarily high frequencies (the product of A’s in 0 adds the corresponding 
wave numbers k)Q 


Theorem 4 Assume that the system ( |25| ) being integrated satisfies the assumptions 
above. Then there exist positive constants h^, C, both independent of tc, such that: 

1. For \h\ < /iQ arbitrary 9 q, the true solution Xg = {yg, 9g), and the 

numerical solution fihixg) are well defined and lie in 12. 


2. The local error at xg satisfies 
fih{xg) - fih{xg) = 


0 

h{a — l)c 


Af-l 

E 

n—l 


h' 


■ E 

ki,....k„Gl 


(-^ki---kn, (^) -^kp-'-kn, (^)) /"kp-'-k^ (^o) 

(46) 


+7^/i(xo), 
where \lZh{xo)\ E C\h\^. 

^For smooth solutions, terms with high frequency must have small amplitude, a fact that may be exploited 
in the derivation of error bounds (Sl.lll. This point will not be studied here. 
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The theorem reduces the estimation of the local error to the estimation of the quan¬ 
tities ylki -k„(^) — ^ki -k„(^)- These are errors arising in the quadrature of scalar 
smooth trigonometric functions and are completely independent of the function f. 

It is assumed hereafter that the integrator is consistent. We first analyze the local 
error in the limit h ^ 0. The condition a = 1 implies that the first term in the right- 
hand side of (46 1 vanishes. Furthermore, from Remark]^ 4lk(/i) — A\^[h) = 0[h) as 
/i —>■ 0 and we conclude that fhix) — 4>h{x) = 0{hf). Note that, for the word with 
0 S as its only letter, Ao{h) — Ao{h) = 0. Moreover, in view of Remark]^ for each 


n < TV — 1, the n-th term in the sum in (461 is actually rather than 0{h'^). 


If, additionally, the underlying univariate quadrature rule is second-order accurate, 
i.e. = 1/2, then A\^{h) — A\^{h) = 0{hf), and the integrator will be second 

order accurate, fhix) — 4>h{x) = 0{hf), provided that hypothesis 3 holds with N > 3. 

The argument may be taken further to translate accuracy properties of the associ¬ 
ated quadrature and cubature rules into accuracy properties of the integrator in the limit 
/i —0. In this way one recovers the order conditions for splitting methods listed in iU 
(cf. ED). We shall not pursue that path: our interest lies in the size of the local error 
when h is not small relative to the periods present in the dynamics, a scenario that we 
discuss next. 

It is well known that for a quadrature rule that is exact for polynomials of degree 

< cr, ^ 

|Ak(/i)-Ak(/i)| <C|k-a;r+ih-+\ 

for a constant C that only depends on the rule. Therefore for the quadrature errors 
to be small it is necessary that \h\ be small with respect to min(l/|k • cli|), where the 
minimum is extended to all k with k • w ^ 0. We reach the unwelcome conclusion that 


the size of the bound in (46 1 depends both on the size of the perturbation / and on oj. 


Example 1 Consider the familiar Strang splitting, r = 2, 

0-1 = 1 / 2 , 02 = ^ 1^1 bi = 1 , &2 = 0 . 


(47) 


The underlying quadrature formula is the second-order accurate midpoint rule. For this 
integrator, for each k such that k • w 0, 

Ak(/i) - Ak(Ti) = exp((l/2)ik • ujh) - —1 (qg) 

IK • ujri 

(for k • a; = 0, Ai^{h) = A]^{h) = 1). An elementary computation leads to the bound 

\Mh)-Ai,{h)\<^\k-uj\^hf 

where the constant 1/24 cannot be improved if the inequality has to hold for arbitrary 
h. □ 


Remark 9 The dependence on uj of the local error is not an artifact introduced by our 
method of analysis. Flere is an example. Consider the forced spring (see the Appendix), 
{d/dt)p = —uj'^q+F, dq/dt = p, where F ^ 0 is a time-independent force and w > 0. 


21 






0.100 
0.010 
0.001 
10 -'' 

10 -= 

0.006 0.01 0.018 0.034 0.064 0.116 0.214 



Figure 1: Energy error at time T = 50 as a function of h (in doubly logarithmic scale) 
for Example]^ The discontinuous straight lines correspond to 0{h^) (left) and 0{h) 
(right). Small circles have been located at points whose abscissa is a value of h that 
leads to a hrst order numerical resonance (Section [573] ). It is apparent that those points 
give rise to local maxima of the error. 


This is the Hamiltonian system with Hamiltonian H = (l/2)p^ + I2)q^ — qF or, 

in action-angle variables 


/ 2a 


1 2a 


1 2a 


H = Lua — \ — sin 9 F = uja — — exp(i0)F + — ejq){—i0)F. 


2i 


2i 


There are two Eourier modes fc = ±1 in the perturbation. 

Choose initial conditions po = 1, go = 0 (with kinetic energy 1/2 and no potential 
energy in the spring). If h/{I/ oj) = 2tt, after one time step, the true solution has 
p{h) = 1 and Strang’s method ( |T^ yields an approximation p{h) = 1 — hF\ therefore 
a bound of the form \p{h) — p{h)\ < C\hY’^^, \h\ < hg, with C and Hq independent 
of ui cannot exist for cr > 0. Note that, after m steps, the error in p will be mh\ 


Example 2 In order to observe the behavior of the Strang’s method ( [47| in problems 
more involved than the scalar example in the last remark, we have integrated the Hamil¬ 
tonian problem with d = 5 degrees of freedom with Hamiltonian function from ll22l 
Chapter XIII.9] 

i=i 

with 

and wi = 1, 012 = 013 = 70, 014 = 70^/2, 015 = 2w2. The system is split by dividing 
the Hamiltonian into its harmonic (quadratic) part, with linear dynamics, and the per¬ 
turbation corresponding to U (see the Appendix). We integrated this problem over the 
(very long) interval 0 < t < 35000 with the initial condition (also taken from 122) 


/ 1 3 7 9 4\ 

V 5’5’5/ ’ 


9o 



3 4 

10u;2 ’ 5 w2 ’ 


llv^ 7 \ 
10w4 ’ ’ 


for which the energy in the harmonic part is 4.225. Eigure [T] gives the error in the 
Hamiltonian at time f = 50 as a function of h. Two different regimes are apparent in 
the figure: 
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1. For h small, the error in the Hamiltonian is very approximately Ch?, as it corre¬ 
sponds to the second order of accuracy of Strang’s splitting. In this regime the 
oscillatory nature of the problem is not relevant and the integrator may be ana¬ 
lyzed by standard techniques, i.e. expansion of the local error in powers ofh and 
transference, using stability, of local error bounds to bounds of the global error. 

2. For h large, the error presents a very irregular behavior. This is due to the highly 
oscillatory character of the solution and, as we shall now describe, may be ana¬ 
lyzed via the word series expansion of the local error. 


5.3 Processing 


It is clear that, as distinct from the global error, the quadrature error in (481 varies 
regularly as h varies. The irregularities in Fig. [T] stem from cancelations, due to the 
oscillations, of local errors in consecutive time steps. For this reason, sharp error esti¬ 
mates in highly oscillatory problems (see e.g. HD, IMI) do not bound the local error 
and then sum the bounds; they rather sum first and bound later, so as to take advantage 
of possible cancelations. We use here an alternative approach that exploits the idea of 
processing that goes back to Butcher G). The presentation here follows f26\ . 

If Xh is a near-identity mapping in and is an integrator, the mapping 


(l^h = Xh Xh 

defines a processed numerical integrator. For to > 1 

= {xh^ °^hO Xh)™ = ° ° Xh; 


(49) 


(50) 


therefore to advance to steps with one may preprocess the initial condition to hnd 
Xh{xo), advance to steps with the original method and then postprocess the numerical 
solution by applying Xh^- Postprocessing is only performed when output is desired, not 
at every time step. In practice, the idea of processing is useful if Xh may be chosen in 
such a way that (f)h is more accurate in some sense than the original (ph'. one then obtains 
extra accuracy at the (hopefully small) price of having to perform the processing (this 
gives rise to Butcher’s notion of effective order Q, JU). Here we use the idea of 
processing as a technique of analysis. We shall process the splitting method ( [40| by 
means of a mapping Xh with an expansion in word series: Xh{x) = WK,(^h){x), K{h) G 
Q (see Remark]^. Then the processed integrator will be expressible as an extended 
word series with coefficients of the form {huj,a{h)) G G- By implication, the local 
error will be of the form (45 i with Akj...k„ {h) — ^ki -k„ (h) in lieu of Aki...k„ (h) — 
^ki -k„(^)- (We have used the obvious notation /i”Aki...k„(^) = Ski -k„ for the 
scaled coefficients of the processed method and will similarly set h"iTki k„(h) = 
trki -k„ (h).) Our policy is to determine the processing, i.e. to determine K{h) G G, in 
such a way that Aki...k„ (h) — Aki...k„ (h) = 0, for all oscillatory words. We then may 
hope that the analysis of the processed integrator would be free from the difficulties 
usually associated with integrators of oscillatory problems. Finally the results on the 
processed integrator obtained in this way will be translated into results for the original 
(ph- 
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5.3.1 First-order numerical resonances 


The conjugation (491 of the mappings may be translated with the help of (34i into the 
equation 

{huj, a(h))it(0, K(h)) = (0, n{h))-k{hu:, a{h)) (51) 

for the coefficients. For words with one letter, ( [5T] l implies, according to the definition 
of ★; 

exp(ik • ujh) K^{h) + A\^{h) = A\^{h) + 

There are two cases to be analyzed. We first look at words (including k = 0) that are 


not oscillatory, i.e. k • w = 0. The value K\^{h) drops from (511 and may be regarded as 
a free parameter. In addition, Ay^{h) = A\^{h) = 1 and therefore A\^{h) — Ay^{h) = 0. 
We then consider oscillatory one-letter words k, k • w 7^ 0, and according to our policy, 
we try to get A\^{h) = Ay{h). This leads to 


K^{h) = 


Ai^jh) - A^jh) 
exp(zk • uih) — 1 ’ 


(52) 


provided that exp(ik • ujh) 7^ 1. If k • w 7^ 0 and exp(ik • ujh) = 1, we say that a 


first-order numerical resonance occurs. When this happens, K\^{h) drops from (511 


and ^k(^) = ^k(^)- As a consequence, in general, Ai^{h) — Ai^{h) will not vanish. 

If, for given h, there is no first-order numerical resonance, then the expansion of 
the local error only contains terms corresponding to words with two or more letters. In 
analogy with Theorem|^(details will not be given), it is then possible to bound the local 
error of the processed integrator by C/i^, with C independent of uj. This in turn will 
lead to a C'h bound for the global error of the processed integrator and, after taking 
into account the pre- and postprocessing to a C"h bound for the global error in the 
method fih being analyzed. The constant C” may be chosen to be independent of h, 
provided that h is bounded away from the resonances; it worsens as h gets closer to a 
numerical resonance in view of (521. This explains the troughs in Fig. [T] 

On the other hand, if for given h there is at least one numerically resonant k S I, 
then, processing is of no help in removing the w-dependent quadrature error of the 
original, unprocessed method. This was only to be expected because at a numerical 
resonance, as shown in Remark]^ the global error in fih may actually be large (cusps 
in Fig. [T). 

Remark 10 By using the operation if to compute the expansion of the m-fold com¬ 
positions (j)]^ and 0™, we find after some simple algebra that if ±1 are numerically 
resonant wavenumbers and exp(ik • ojh) 7^ 1 for k 7^ 0, ±1, then the error over m 
steps has an expansion 


C(^o)-CM = 

mh (^Ai{h)fi{xo) + A_i(/i)/_i 

exp(ik • uimh) — 1 
exp(*k • ujh) — 1 


E 


kGl\{0,±l} 


(Ak(/i) - Ak(/i))/k(a:o) 
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Thus the mh growth as m increases with fixed h we already encountered in Remark|^ 
holds for general integrators and general differential equations. 


5.3.2 Higher-order resonances 


Assuming that h does not satisfy any first-order numerical resonance, one may go a step 
further and look at words with two letters kl; these are oscillatory if (k -f 1) • w 7^ 0. 
Now •HD implies 

exp(z(k -I- 1) • ojh) K]f^i{h) + Ai^{h) exp(jl • ujh)K\{h) + Aki(fi) 

= ^ki(^) + Kkih)Ai{h) + Kki{h). 


Whenever exp(i(k + 1) • ojh) = 1 (second order numerical resonance), the value of 
K\a{h) cannot be chosen to ensure that Aki(/i) = Aki(fi). A similar consideration 
applies for nonoscillatory words with n letters when exp(i(ki -f • • • -f k„) • ujh) — 1. 
When no numerical resonance takes place, the equation ( [ST] ) may be used to find 
values Ku]{h), w G W such that, on the one hand, define an element K{h) that be¬ 
longs to Q, (i.e. the shuffle relations hold) and, on the other, ensure that Aki k„ (h) — 
^ki -k„ (h) — 0, for all oscillatory words. This will be proved in Remark 13 below by 
using modified systems. 


Remark 11 Since pre- and postprocessing introduce in any case 0{h) errors, the pro¬ 
cessing technique used here yields 0{h) bounds for the global error of (jih even for 
values of h where there is no first-order or second-order numerical resonances. Fig. 
shows that, for this simulation, 0{h?) global error bounds cannot exist if h is large 
relative to the periods present in the dynamics. 


5.4 Modified equations and modified Hamiltonians 

Modified equations ll20l . 191 , l35l . l22l provide a useful means to describe the be¬ 
haviour of numerical integrators. 


5.4.1 Modified system using one letter words 


We look for a (one-letter word) modified system 


= W 
dt 


(ui,P(h)) 




(53) 


where /3 S g, /3uj = 0 for w G W„, n > 1 and the coefficients /3k(h), k € X are chosen 
in such a way that, for words with one letter, the extended word series expansion of 


the fi-flow of (|53| matches the corresponding expansion for the integrator (f>h- We 
recall that the system being solved is also of the form (53 i with the coefficients /3 given 
in I 


By integrating the system ( [53] l by the procedure outlined in Section|^and imposing 
that its flow matches (j)h to the desired order, we find the condition 

exp(ik • ujh) — 1 


ik • ujh 


Mh) = Ak(fi) 


(54) 
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(it is understood that the fraction takes the value 1 if k • a; = 0). For k = 0 or for 
any one letter word that is not oscillatory, this implies /3\^{h) = 1. For an oscillatory 
one-letter word k ^ 0, k € I, if ft, is such that ik • ujh — 211 j for some integer j ^ 0 


(first order numerical resonance), then the fraction in (54i vanishes and the equation 
for /3k{h) will in general not be solvable. Thus first-order numerical resonances are 
obstructions to the construction of the modified system. (In fact the counterexample 
in Remark Improves that modified systems of the form envisaged here do not exist at 
numerical resonances.) When ft is bounded away from resonances, — ph = 0{h?) 
with the implied constant independent of w (as in Theorem]^. 


Example 3 For Strang’s method, if k is oscillatory and there is not a numerical reso¬ 


nance, (541 yields the value 

Mh) = 


k • (jjh 


2 sin(k • ujh/2) 


□ 


Remark 12 Assume that, for given ft, the modified system above has been found. We 
may then try to find a change of variables x = so that in the new vari¬ 

ables the modified vector field matches the field W for words with one letter. 
According to Section]^ we have to impose that /3 * K(ft) + and K{h)-kl3{h) 

coincide for words with one letter. This leads to 

i(k • w) Kk{h) -I- 1 = (3k{h). 

If k is not oscillatory, Kk{h) is free because, as noted above, /3k(^) = 1- For k 


oscillatory, Kk(ft) is uniquely determined. By using (54i, a little algebra shows that the 


value of Kk(ft) found in this way is the same we obtained in (52i. Thus the change of 


variables WK(h) we used for processing may be seen as determined by the requirement 
that, in the new variables and for nonoscillatory one-letter words, the modified vector 
field of the unprocessed integrator reproduces the vector field being integrated. 


5.4.2 Other modified systems 

More precise modified systems may be constructed by successively adding to the mod¬ 
ified vector field contributions from words of 2, 3, ... letters. For the n-th of these 
modified systems, the modified vector field has Pui = 0 for words with more than n 
letters and we impose that, for words with n or fewer letters, the extended word series 
expansion of the ft-flow matches the corTesponding expansion for the integrator 
(ph- 

For two-letter words, proceeding as in the case of one-letter word modified systems, 
we obtain the condition 


exp{i(k -f 1) • ujh) — 1 
i(k -f 1) • (jjh 


/3ki(ft) + ft4lki(ft)/3k(/r)A(M = hA^,{h). 


In general, the resulting equation is of the form 

exp(i(ki + • • • 4- k„) • wft) — 1 


*(ki -f • • • k„) • wft 


/3ki...k„(/i) - ft"-'Ak,...k„(/r) = 0(ft”-'), 
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where the right-hand side depends polynomially on the coefficients f3w{h) of words 
w with less than n letters. Such an equation can be solved for /3ki -k„(^) provided 
that there is no numerical resonance, (ki -I- • • • -f k^) • ooh = 27 rj, r < n, j ^ 0. In 
the limit where the length of the words increases indefinitely one obtains, if there is 
no numerical resonance of any order, a modified system whose formal fi-flow exactly 
reproduces the expansion of (ph- 

As explained in Section]^ the modified systems found in this way are Hamilto¬ 
nian whenever the system being integrated is Hamiltonian. Furthermore the modified 
Hamiltonian functions are easily expressible in terms of brackets. 


Example 4 For the linear forced oscillator in Remark each word basis functions 
associated with words with two or more letters vanishes. In this case the modified 
systems above with n > 1, coincide with the modified system using only one-letter 
words and the latter is exact, ie = <ph- The (exact) modified Hamiltonian is 

and therefore in the particular case of Strang’s method we have 


w 


Loh 


2” +¥« - 1 2An{u,hn) 


F. 


For nonresonant h, the effect of using the splitting method is to alter the value of the 
applied force. Unless \ujh\ ^ 1 the misrepresentation of the force introduced by the 
discretisation will be large. We emphasize that, as distinct from the situation when 
using conventional modified equations based on series of powers of h, the analysis 
here does not require h to be small. □ 


Example 5 For the problem in Example we have measured the variation as t in¬ 
creases of the true energy H of the numerical solution and of the corresponding varia¬ 
tions of the energies in the modified one-letter-word Hamiltonian and two-letter-word 
Hamiltonians. We used h = 0.7974 as this, which is more than 12 times larger than 
the period of the fastest oscillator, avoids first and second order resonances. The results 
given in Figs. |2]^ clearly bear out how the one-word-letter modified system matches 
the numerical solution much better than the system being integrated, but not as well as 
the two-word-letter modified system. 


Remark 13 The idea in Remark 12 may be extended. Assume that there is no nu¬ 
merical resonance of any order so that it is possible to find a modified equation whose 
formal flow exactly reproduces the expansion of the integrator. By using Theorem 
we may bring the modified system to a normal form where the contribution of all oscil¬ 
latory words have disappeared. This implies that a processor has been found such that 
the expansion of the local error of the processed method does not contain contributions 
of oscillatory words. 
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Figure 2; Variation of the true Hamiltonian (energy) evaluated at the numerical solution 
as a function of t for Example|^(/i = 0.7974). 
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Figure 3; Variation of the one-letter-word Hamiltonian evaluated at the numerical so¬ 
lution as a function of t for Example|^(/i = 0.7974). The vertical scale is 100 times 
larger than in Eig.[^ 
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Eigure 4; Variation of the two-letter-word Hamiltonian evaluated at the numerical so¬ 
lution as a function of t for Example|^(/i = 0.7974). 
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We close this section with an observation. As noted before, numerical resonances 
((ki + • • • + kr) • ojh = 27rj, j ^ 0 ) obstruct the construction of modibed systems; 
nonoscillatory words ((ki + • • • + k^) • w = 0) cause no trouble in that connection. On 
the other hand, the nonoscillatory character of a word is an obstruction to its elimination 
by changing variables. Processing, that as we have just seen is equivalent to hnding 
modibed problems and then changing variables, is hampered by both numerical res¬ 
onances and non-oscillatory terms. See in this connection ( [52] i, whose denominator 
vanishes both at numerical resonances and for nonoscillatory words. 

6 Technical results 

This section is devoted to more technical material. 


6.1 Algebraic results 

6.1.1 Differential operators 

In ([^ we associated with each vector held in Q a brst-order linear differential 
operator Ea- With each word w = ai ■ ■ ■ an, n > 0 we now associate the n-th order 
(linear) differential operator obtained by composition; 


Ea-,...a„g{x) = Eai ■ Ea^ ■ ■ ■ Ea^g{x)\ 

Efh is debned as the identity operator. Finally, with each 7 S we associate the 
formal series of linear differential operators 


D,y = 






Two non-empty words w = ai ■ ■ ■ am, w' = a'^ - ■ ■ a'n may be concatenated El 
to give rise to the word ww' = oi • • • amO-'i ■ ■ • In addition %w = w% = w for each 
w G W. Clearly, concatenation of words corresponds to composition of the associated 
operators: Emw' = E^] ■ Em'- 

Each word ui • • • a„ may be deconcatenated in n + 1 different ways: 0(ai • • • a„), 
(ai)(a 2 • • • a„),..., (oi • • • a„)0; these feature in the debnition (11 1 of the convolution 
product. This observation leads to the following rule for the composition of two series 
of operators; 


■ Dy = D. 




(55) 


6.1.2 The shuffle algebra 

The product lu may be extended in a bilinear way from words to linear combinations 
of words, i.e. if pj/ are scalars; 

(XI ^ (X = X ^ ■ 

3 3 ' 3 , 3 ' 
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When endowed with this operation, the vector space C(y4) of all such linear combi¬ 
nations is a unital, commutative, associative algebra, the shuffle algebra, denoted by 
sh(A) (see OTl . ifTSl . ESll '). Note that sh(A) is graded by the number of letters of the 
words. 

Deconcatenation defines a coproduct and turns sh(A) into a (commutative, con¬ 
nected, graded) Hopf algebra 0. It is well known that the dual vector space of a Hopf 
algebra is automatically endowed with a product operation. Here the dual of sh(A) 
may be identified in a natural way with by associating with each linear form i on 
sh(A) the family of coefficients 7^, = (.{w), w G W. After this identification, the 
product in the dual of sh(A) coincides with the convolution product * defined in (11 1 . 
The sets Q and g in Section]^ are then respectively the group of characters and the Lie 
algebra of infinitesimal characters of the Hopf algebra sh(A); well known results on 
Hopf algebras show that exp^ in ( [T6| maps g onto Q and has an inverse given by log,^ 
in ( [T7 | i, see e.g. Il28l . ifTSl . 


6.1.3 The actions of G and g on word series 

As shown e.g. in lfT4l . there is a narrow connection between the word basis functions 
fw{x) and the operators E^, w GW: 

fw{x) — E^x, 

(in the right-hand side, with an abuse of notation, x denotes the identity function that 
maps each Z3-vector into itself). As a consequence we have the following correspon¬ 
dence between word series and series of operators 

W^{x) = D^x. (56) 

The use of series of operators is common in control theory and dynamical systems; 
word series, being series of functions, provide a more convenient way to study numer¬ 
ical integrators. 

The operators Ea, a G A, are derivations: Ea{gh) = {Eag)h -f g{Eah) for each 
pair of scalar functions g, h. Iteration yields: 

Eabigh) = {Eabg){E(/)h) + {Ebg){Eah) + {Eag){Ebh) + {E(/)g){Eabh), 

Eabcigh) = {Eabcg){E(iih) + {Ebcg){Eah) + {Eacg){Ebh) + {Ecg){Eabh) + 

{Eabg){Ech) + {Ebg){Each) + {Eag){Ebch) + {E$g){Eabch), 

etc. Note that, in the hrst of these identities, the pairs of words {ab, 0), (a, b), {b, a), 
{fti,ab) that feature in the right-hand side are precisely those whose shuffle product 
gives rise to the word ab that appears in the left-hand side. A similar observation may be 
made in the second identity. In general, if ly € Wm, w' G Wn and wluw' = xoj, 
then the Wj G Wm+n are precisely those words for which {Eyjg){Eyjih) is one of the 
2 m+n Qf j-jjg expansion of This result may be used in combination 

with the shuffle relations ( [T^ to prove (see e.g. ifTSlI . Theorem 2 ) that, for j G G, 

Ejigh) = D^{g) D^{h). 

^Algebraically, the action of the operators Ew on products gh defines a coproduct (6); the shuffle product 
is obtained from this coproduct by duality (m Section 1.5]. 
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By considering the coordinate mappings g{x) = x^, h{x) = 
that 


and (56 1 we conclude 


x\W^{x)) x\W^{x)) = D^{x^)D^{x^) = D^{x\^) 


and then linearity shows that, for each polynomial mapping P, P{Wj{x)) 
It follows that 


g{Wj{x)) = Djg{x), l&Q] 


D^P{x). 

(57) 


for any (scalar or vector valued) smooth mapping g. Thus D^g{x) provides the formal 
expansion of the composition g{W.^{x)) provided that the coefficients 7 belong to the 
group Q. 

The proof of the formula ( [T3] l, that defines an action of the group Q on the vector 
space of all word series, is now easy: 


Ws{W^{x)) = D^Ws(x) = {Dy ■ Ds)x = 


we have successively used ( |57] i, ( |56l l, ( |55l l and once more 1 

In (18 1 , the expression {dxWs{x))Wij is the result of applying to the word se¬ 
ries Ws(^, the first-order differential operator associated with the formal vector held 
W/ 3 (x). The formula then reveals that the action of the algebra g on word series corre¬ 
sponds to the operation /3 * <5. 


6.1.4 Linear differential equations 

In SectionSit was proved that the initial value problem (j2T|l has a unique solution with 
a{t) G for each t. We show here that in fact a{t) G Q. We use the following 
auxiliary result: 


Lemma 1 Assume that ry G is such that, for some positive integer n, and for each 
G yVi, w' G yVm, l + m < n, the shuffle relation ( 72 I hold. Then, for (3 G Q, 


w 


w G Wi, w' G yVm, I + m < n, with w luw' = 


Wi 


^ ^ — Vw ijl I3fw' ij] ^ fffwXjw' • 

3 

Proof: Since j3 G q, there exists a curve 7 (f) in Q such that ( |T5| ) holds. The hypothesis 
of the lemma then allows us to write: 




The result is obtained by applying d/dt\^^Q to both sides of this equality. □ 


Now consider the solution a{t) G of (21 1 . We shall prove by induction on n 
that for each t 

^ ^ c^{f)wj — (58) 

3 

for w G Wi, w' G Vdm, I + m < n, with wluw' = J2j '^3- 
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This trivially holds for n = 0, since a(f )0 = 1 for all t. Assume that (58 1 is 
satished for some n > 0 , and choose w G Wi, w' G Wm, I + m < n + 1 , with 
WLUw' = From (211 we hnd 


d 

dt 


* /3{t))^. - (a(f) */?(f))u,a(f)u,/ - a{t)w{a(t) * |3{t))^ 


and the lemma implies that the right-hand side of this equality vanishes. Since (58i 
holds at f = 0 , it does so for each value of t. 


6.2 Proof of Theorem [T] 


We simplify the system ( [3^ by performing a sequence of changes of variables with 
coefficients ...in ^ so that the change dehned by simplihes the 

coefficients of the vector field associated with words with n letters and leaves unaltered 
the coefficients associated with shorter words. The element is sought in the form 
exp^(A*^[”F) where g g and aIS"^^ = 0 if u; G W \ yV„. If and /Sl"! 

are respectively the vector helds before and after the n-th change of variables, and 
u; = ki • • • k„, the equation (38 1 implies, after taking into account that vanishes 


for nonempty words with less than n letters 


i((ki H- kn) ■ uj)kw =/3^w 


If (ki + • • • k„) • w 7 ^ 0 we may choose aII^^ to enforce (3^^ = 0. In other case, we 
set = /3w~^^ and aII^^ = 0. The element A^"! constructed in this way belongs to 
g because the required shuffle relations hold (if shuffling two words leads to resonant 
words all the coefficients in A^"! vanish; in the nonresonant case the coefficients AI"! 
are proportional to the corresponding coefficients in which satisfy the shuffle 

relations). In turn G g because 

pin] ^ ^[n] ^ ^[n-1] ^ ; 


both terms of the right-hand side are in g (the second is the value at f = 0 of 

and, as noted above G G). 


6.3 Proof of Theorem m 

For (y, 6) G the function / in ( p5| ) is bounded and Lipschitz continuous. Therefore, 
for |f| small {y{t),9(t)) = 4>tiyoi9o) is well defined and |y(t) — yo| < Ci\t\, where 
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Cl is a bound for |/|. A simple contradiction argument shows that \y{h) — t/ol < for 
\h\ < R/Ci. 

To deal now with the numerical solution, dehne the intermediate points (stages), 
J = 1 , ■ ■ ■ C, 


If \bjh\ < R/{Ci/r), the iteration of the argument used above ensures that \yj — 
yj-i\ < R/r, j = 1,... ,r and then the triangle inequality implies that (j)h{yo, ^o) G 

n. 

_fjV) _ (N) 

We shall use the notations refer to the result of 

suppressing all terms corresponding to words with N or more letters of the extended 
word series with coefficients a{t), a{h) respectively (of course the alphabet is now X 
rather than Z‘^). In addition we set 


7^ 


(T)/ 


h {x) = - W 


TjyW ! 1 

{huj,a.{h)) \^)■> 


7^ 




{hiij,a{h)) 


(x) 


(the superscripts T and S mean ‘true’ and ‘splitting’). Our task is to bound TZh{xo) = 
(a;o) —(xq). For (xq), by stopping the iterative procedure (see e.g. lIT^ l 
that leads to ([28]l we find the following representation: 





dtN exp('ikAr • wfAf) • • • 



dti exp(iki • a;fi)/ki...k„(^ii(a;o))- 


We know that (j)ti{xo) G 0 for |/i| < hg, and therefore we may guarantee that 
\'RP\xg)\ < C\h\^, with C depending only on X and bounds for the derivatives 

of the Fourier coefficients. _ _ 

For TZP{xg) we use a similar device. The key point is that (cf. (281-(29 1 ) 


(khiyo.dg) = (0, h(ai 


■ar)uj) + {y{h),r]{h)), 


where {y{t),r]{t)) is the solution of 


^ \y 

dt rj 


Ak(f) h{y,v), 

kGZ'i 


with piece-wise constant functions defined by 

Ak(f) = xbj exp (zk • uj{ai + • • • + aj)h), {j — l)h/r < t < jh/r, 1 < j < r. 


This differential system associated with (j)h is very similar to the system (271 associated 
with (ph, the difference being that for the former the complex exponentials are frozen at 
the times (oi + • • • + aj)h. After this observation the residual for ph is bounded with 
the technique used for the true ph- 
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Appendix: examples of perturbed integrable problems 

Any system 

= Mw + F{w), (59) 

where M is a skew-symmetric D x D constant matrix, may be brought by a linear 
change of variables to the form: 

id/dt)z^ = r(z,P,Q), l<j<D-2d, (60) 

(d/dt)P^ = + g\z,P,Q), 1 < f < d, 

{d/dt)Q^ = P^ + h\z,P,Q), l<i<d 


Here d is the number of nonzero eigenvalue pairs izwr, ui > 0, of M. In the unper¬ 
turbed case, P = 0, = 0, = 0, the system consists of d uncoupled harmonic 

oscillators with frequencies together with D — 2d trivial equations {d/dt)z^ = 0. 
The introduction of the variables a^, 9^ such that 


= V2ujp cos6i^ 



l<£<d, 


(61) 
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takes now the system to the format (251 with y = (z^, 


^D-2d 


, a). The system (59 1 


or (60 1 is a natural candidate to integration by splitting methods based on separating the 
linear part (that may be integrated in closed form) from the perturbation. The later may 
perhaps be treated by means of a numerical integrator with a very hne time step. In 
favorable instances, the perturbation may be integrated analytically in closed form; this 
is the situation in the following particular case of ( [60| , commonly found in mechanics 
{D is even and z = (p, q)), 


{d/df)ff = 

f{q,Q), l<j<D/2-d, 

(62) 

{d/dt)q^ = 

ff, 


{d/dt)P^ = 

-f/(g,Q), l<i<d, 


{d/df)Q^ = 

P^. 



Under the dynamics of the perturbation, p and P remain constant and q and Q grow 
linearly with t. 

The system (jh^ is Hamiltonian if the forces f^, derive from a potential. When 
that happens, the introduction of the canonical {dP^ A dQ^ = da^ A dO^) action/angle 
variables in ( |6T] l preserves the Hamiltonian character of the equations of motion and 
even the value of the Hamiltonian function. 

So far the unperturbed problem has been linear, but nonlinear cases may also be 
treated. Typically, integrable nonlinear problems may be brought to the form ( [4l] i with 
u) = oj{y)-, a device commonly used e.g. in dynamical astronomy consists in hxing a 
relevant value yo of y, decomposing uj{y) = uj{yQ) + A(j/) and seeing (0, A(y)) as 
part of the perturbation. 

There are many instances of perturbations of nonlinear integrable systems, after 
the introduction of suitable action/angle variables take the form ( [25] l. A well-known 
example is provided by perturbations of the Keplerian motion of a celestial body. 


Remark 14 For ( [62] i, as we just noticed, the system corresponding to the perturbation 
may be solved in closed form in the variables (p, q, P, Q). On the other hand, the anal¬ 
ysis in Section 1^ operated with a different set of variables x = {y,0) = {p,q,a,6). 
This causes no difficulty: it is standard practice when using splitting methods that the 
different split systems are integrated employing different sets of dependent variables. 
In partial differential equations, parts corresponding to linear, constant-coefficient dif¬ 
ferential operators are typically integrated in Fourier space and nonlinearities in phys¬ 
ical space. Splitting methods are based on true solution flows, which of course com¬ 
mute with changes of variables. The situation is very different for, say, Runge-Kutta 
schemes, where (except for affine changes) changing variables does not commute with 
the application of the numerical method, and the performance of the integrator very 
much depends on the choice of dependent variables. For instance, any consistent 
Runge-Kutta method integrates exactly the unperturbed problem when written as in 
(41 1 but incurs in errors when dealing with the unperturbed version of (601 (/* = 0, 
= 0 , = 0 ). 
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